Below is an example of 20 hypothetical 1ha restoration plots (10,000m2, 100m*100m) distributed at random at three reefs in the Townsville region (Davies Reef, Big Broadhurst Reef, Little Broadhurst Reef). The code takes the Allen Coral atlas geomorphology data (see here for more details), filters all reef polygons larger than 1ha in size, randomly selects 20 polygons across the three reefs, and constructs a 100*100m plot around the centroid of each polygon:

# data packages
library(tidyverse)
library(units)
library(foreach)

# spatial packages
library(sf)

# mapping packages
library(tmap)
library(leaflet)

# create function to get make polygons for plots
set_restoration_plot_centres <- function(input = NULL, width = NULL, length = NULL) {

  # Calculate the coordinates of the rectangular polygon
  x <- sf::st_coordinates(input)[1, 1]
  y <- sf::st_coordinates(input)[1, 2]

  # set parameters
  x_min <- x - (width / 2)
  x_max <- x + (width / 2)
  y_min <- y - (length / 2)
  y_max <- y + (length / 2)

  polygon <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
    sf::st_sfc(crs = 20353)

  return(polygon)
}



habitats <- c("Plateau", "Back Reef Slope", "Reef Slope", "Sheltered Reef Slope", "Inner Reef Flat", "Outer Reef Flat", "Reef Crest")

davies_geomorphic <- st_read("/Users/rof011/coraldynamics/data/Davies-Broadhurst-20230906193843/Geomorphic-Map/geomorphic.geojson") %>%
   sf::st_transform(crs = sf::st_crs("EPSG:20353")) %>%
   dplyr::group_by(class) %>%
   dplyr::mutate(habitat_id = paste0(
   gsub(" ", "_", class), "_",
   sprintf(paste0("%0", ceiling(log10(max(1:length(class)))), "d"), 1:length(class)))) %>%
   sf::st_make_valid() %>% 
#  mutate(habitat_id=as.factor(habitat_id)) %>%
  dplyr::group_by(habitat_id, class) %>% 
  summarise(geometry = st_union(geometry)) %>% 
  mutate(area = round(st_area(geometry))) %>%
  filter(class %in% habitats) %>%
  drop_na(class) 
## Reading layer `Davies Broadhurst' from data source `/Users/rof011/coraldynamics/data/Davies-Broadhurst-20230906193843/Geomorphic-Map/geomorphic.geojson' using driver `GeoJSON'
## Simple feature collection with 2291 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 147.6256 ymin: -18.98455 xmax: 147.8017 ymax: -18.80138
## Geodetic CRS:  WGS 84
davies_plots_centroids <- davies_geomorphic %>%
  filter(area > set_units(10000,"m^2")) %>%
  filter(class %in% c("Reef Slope", "Sheltered Reef Slope", "Back Reef Slope")) %>% 
  drop_na(class) %>%
  st_centroid() %>% 
  st_cast("POINT")

library(foreach)

davies_plot_outputs <- foreach(i=1:nrow(davies_plots_centroids), .combine="c") %do% {
  tmp <- set_restoration_plot_centres(davies_plots_centroids$geometry[i], width=100, length=100)
  tmp
}

Map restoration plots

Select layers using the layercontrol on the left, use [ ] for full-screen viewing.

library(tmap)
set.seed(1)
davies_plot_outputs <- davies_plot_outputs |> st_as_sf() |> mutate(id=paste0("Plot_ID_",seq(1,n(),1))) 

davies_plot_outputs2 <- davies_plot_outputs[sample(nrow(davies_plot_outputs), 20), ]


davies_plot_outputs <- davies_plot_outputs |> st_as_sf() |> mutate(id=paste0("Plot_ID_",seq(1,n(),1))) 

davies_plot_outputs2 <- davies_plot_outputs[sample(nrow(davies_plot_outputs), 20), ] # subset to 20

# visualise with thematic maps
map <- tm_view() +
tm_tiles("Esri.WorldImagery", group = "<b> [Seascape]</b> satellite map", alpha = 0.5) +
  tm_shape(davies_geomorphic |> mutate(class=as.factor(class)), id="area", name = "<b> [Seascape]</b> habitats") +
  tm_borders(col = "black", lwd = 0.2) +
  tm_fill("class", name = "Benthic habitats", id="area", palette = c("Plateau" = "lightskyblue1", "Back Reef Slope" = "darkcyan",
                                                                    "Reef Slope" = "darkseagreen4", "Sheltered Reef Slope" = "darkslategrey",
                                                                    "Inner Reef Flat" = "darkgoldenrod4", "Outer Reef Flat" = "darkgoldenrod2",
                                                                    "Reef Crest" = "coral3"), alpha = 0.6) +
  tm_shape(davies_plot_outputs2, id="plot_id", name = "<b> [Restoration]</b> plots") +
    tm_fill(col="darkred", alpha=0.6) +
    tm_borders(col="black")

  map |> tmap::tmap_leaflet() |>
    leaflet::addProviderTiles('Esri.WorldImagery', group = "<b> [Seascape]</b> satellite map", options=leaflet::providerTileOptions(maxNativeZoom=18,maxZoom=100)) |>
    leaflet::addProviderTiles('Esri.WorldTopoMap',  group = "<b> [Seascape]</b> base map", options=leaflet::providerTileOptions(maxNativeZoom=19,maxZoom=100)) |>
    leaflet::addLayersControl(position="topleft", overlayGroups=c(
      "<b> [Seascape]</b> satellite map", "<b> [Seascape]</b> habitats", "<b> [Restoration]</b> plots"),
                              options=leaflet::layersControlOptions(collapsed = FALSE)) |>
    leaflet.extras::addFullscreenControl(position = "topleft", pseudoFullscreen = FALSE)

Restoration area

What does 20ha cover across the three reefs? Using the spatial data we can calculate restoration area as a proportion of reef habitat (combined reefs: Davies, Big Broadhurst, Little Broadhurst) assuming 20ha is located across i) Total coral area, ii) Total reef slope area, iii) Total reef flat area:

Habitat Area Proportion
Total coral area 35139452 0.6%
Total reef slope area 21021120 1%
Total reef flat area 14118332 1.4%